Librerías requeridas: dplyr, flextable, tidyr, ggplot2, DoE.base2 Impacto de la eliminación de réplicas y puntos centrales en diseños factoriales fraccionales \(2^{k-p}\)
Nota
El presente documento corresponde a un apartado del trabajo titulado “Guía para la implementación de diseños experimentales en la industria farmacéutica”, presentado por Nathalia Cortés Duque como requisito para optar al título de Magíster en Ciencias - Estadística de la Universidad Nacional de Colombia.
3 Metodología
3.1 Evaluación de la significancia de los efectos
El presente estudio tiene como propósito analizar el impacto de eliminar réplicas o los puntos centrales en un diseño factorial fraccional \(2^{k-p}\) sobre la significancia estadística de los efectos principales e interacciones, así como sobre el error cuadrático medio (ECM). Para este fin se realizaron 10 simulaciones independientes, y en cada simulación se generó un caso particular de efectos activos y su correspondiente respuesta (Ver Metodología aquí).
A partir del diseño factorial completo \(2^k\) con dos réplicas y cinco puntos centrales (dis0_cod) con respuesta asignada, se eliminaron las filas correspondientes para obtener un diseño fraccional \(2^{k-p}\) que conserva las dos réplicas y los cinco puntos centrales (dis1_cod). Adicionalmente, se generaron tres variantes eliminando aleatoriamente filas de réplicas o puntos centrales: un diseño con réplicas pero sin puntos centrales (dis2_cod), un diseño sin réplicas con cinco puntos centrales (dis3_cod) y un diseño sin réplicas con únicamente dos puntos centrales (dis4_cod).
Para cada diseño se ajusta el modelo aditivo y se calcula el ANOVA correspondiente. Los resultados de significancia de los efectos se almacenan en las tablas TAB1, TAB2, TAB3 y TAB4, donde el valor 1 indica significancia estadística (\(p < 0.05\)) y 0 indica no significancia. Este procedimiento se implementa evaluando, en cada fila del ANOVA, si el valor \(p\) es menor o igual a 0.05, transformando el resultado lógico en un valor entero.
A partir de estas tablas se calcula el porcentaje de concordancia entre los efectos reales y los detectados por cada modelo, la potencia (proporción de efectos verdaderamente significativos correctamente identificados) y el error tipo I (proporción de efectos nulos detectados como significativos).
3.2 Evaluación del ECM con medidas repetidas
Con relación a la evaluación del ECM, se utiliza la metodología de medidas repetidas, en la cual a una misma unidad experimental se le aplican sucesivamente varios tratamientos, generando valores repetidos de la respuesta sobre el mismo individuo (Dı́az Monroy y Mario Alfonso (2012)).
En este estudio, a partir de los diseños simulados en la sección anterior, cada fila corresponde a la simulación del diseño experimental de referencia (dis1_cod), al que se le generan variantes al eliminar réplicas o puntos centrales (dis2_cod, dis3_cod, dis4_cod). Para cada diseño se obtiene como variable respuesta el Error Cuadrático Medio (ECM), calculado mediante la siguiente fórmula:
\[ \text{ECM} = \frac{\text{deviance}(modelo)}{\text{grados de libertad del residual}(modelo)}. \]
Estos resultados fueron compilados en la matriz ECM_tab, que corresponde a la matriz \(X = (x_{ij})\), donde \(x_{ij}\) representa la respuesta del \(i\)-ésimo sujeto al \(j\)-ésimo tratamiento. Incialmente, se realiza un análisis descriptivo de los ECM (Ver Figura 5.1). Luego, se construye la matriz de contrastes \(C\), donde cada fila representa una hipótesis a evaluar. En este caso, se plantea la utilización de una matriz de contrastes \(C_1\) cuyo objetivo es contrastar simultáneamente todas las medias de los tratamientos. Con ello, se busca determinar si existe al menos un tratamiento que difiera significativamente de los demás, lo que constituye el análogo a la prueba ANOVA clásica, pero en un contexto donde los tratamientos son dependientes.
\[ C_1 = \begin{bmatrix} -1 & -1 & 0 & 0 \\ 0 & -1 & 1 & 0 \\ 0 & 0 & -1 & 1 \end{bmatrix} \]
Asimismo, es posible considerar una matriz de contrastes generales de interés \(C_2\), tales como:
\[ C_2 = \begin{bmatrix} -1 & -1 & 1 & 1 \\ -1 & 1 & 1 & 1 \\ 1 & -1 & 0 & 0 \end{bmatrix} \]
Donde los contrastes a evaluar son:
Contraste 1: Contraste que relaciona los diseños con réplica con respecto a los diseños sin réplica \(C_1 =\) (-1, -1, 1, 1)
Contraste 2: Contraste que relaciona el diseño con réplica y 5 puntos centrales frente al conjunto de variantes derivadas del mismo \(C_2 =\) (-1, 1, 1, 1)
Contraste 3: Contraste que relaciona el impacto de la réplica en diseños con la misma cantidad de puntos centrales \(C_3 =\) (1, -1, 0, 0)
La matriz \(C\) se utiliza para transformar los datos: \(Y = CX\). De esta forma, el problema se reduce a \((p-1)\) dimensiones y se trabaja con contrastes independientes, evitando dependencia entre medias y permitiendo aplicar la prueba \(T^2\) de Hotelling.
El cálculo implementado en el código sigue la formulación teórica:
\[ T^2 = n(C \overline{X})' (CSC')^{-1} (C \overline{X}), \]
donde \(n\) es el número de simulaciones y \(S\) la matriz de covarianzas estimada mediante
\[ S = \frac{1}{n-1}X'(I - \tfrac{1}{n}\mathbf{1}\mathbf{1}')X. \]
La decisión es rechazar \(H_0\) si, para un nivel de significancia \(\alpha = 0.05\),
\[ F = \frac{n - p + 1}{(n-1)(p-1)} T^2 > F_{(\alpha, p-1, n-p+1)} \]
4 Código
4.1 Diseño \(2^{3-1}\)
4.1.1 Significancia de los efectos
Código
# Matrices
TAB1_23 <- matrix(0, nrow = n, ncol = 3)
TAB2_23 <- matrix(0, nrow = n, ncol = 3)
TAB3_23 <- matrix(0, nrow = n, ncol = 3)
TAB4_23 <- matrix(0, nrow = n, ncol = 3)
efectos_reales_23 <- matrix(0, nrow = n, ncol = 7)
names_23 <- c("Ap","Bp","Cp") # Los que serán estimados por el modelo
colnames(TAB1_23) <- names_23
colnames(TAB2_23) <- names_23
colnames(TAB3_23) <- names_23
colnames(TAB4_23) <- names_23
colnames(efectos_reales_23) <- c("A","B","C","AB","AC","BC","ABC") # Todos
# Tabla para completar los errores
ECM_tab_23 <- data.frame(
dis1 = numeric(n),
dis2 = numeric(n),
dis3 = numeric(n),
dis4 = numeric(n))
# Bucle principal
set.seed(3)
for(i in seq_len(n)){
## Efectos a evaluar en el este diseño
todos_efectos_23 <- c("A","B","C","AB","AC","BC","ABC") #Todos
# Selección de un número aleatorio entre 1 y 7
# dado que existen 7 efectos posibles a evaluar en este diseño.
n_ef <- sample(0:length(todos_efectos_23), 1)
# Elección aleatoria _n_ef_ efectos que serán los efectos significativos
efectos_activos <- sample(todos_efectos_23, n_ef)
efectos_reales_23_vec <- setNames(rep(0, length(todos_efectos_23)), todos_efectos_23)
# Asignación de un valor aleatorio a los efectos elegidos como significativos:
efectos_reales_23_vec[efectos_activos] <- sample(c(runif(n_ef, -7, -0.5),
runif(n_ef, 0.5, 7)), n_ef)
# Guardar efectos reales como 1/0
efectos_reales_23[i,] <- c(
as.integer(efectos_reales_23_vec["A"] != 0),
as.integer(efectos_reales_23_vec["B"] != 0),
as.integer(efectos_reales_23_vec["C"] != 0),
as.integer(efectos_reales_23_vec["AB"] != 0),
as.integer(efectos_reales_23_vec["AC"] != 0),
as.integer(efectos_reales_23_vec["BC"] != 0),
as.integer(efectos_reales_23_vec["ABC"] != 0))
# Generar respuesta
## Diseño factorial completo con réplicas y con 5 puntos centrales
dis0_cod <- data.frame(
A = c(rep(c(-1,1),4), rep(c(-1,1),4), rep(0,5)),
B = c(rep(c(rep(-1,2), rep(1,2)), 4), rep(0,5)),
C = c(rep(-1,4), rep(1,4), rep(-1,4), rep(1,4), rep(0,5)))
n_corridas <- nrow(dis0_cod)
## Generación del error
error <- rnorm(n_corridas, 0, 2) #error ~ (mu = 0, sigma^2 = 2^2)
## Generación de la respuesta
yrta <- 50 +
efectos_reales_23_vec["A"] * dis0_cod$A +
efectos_reales_23_vec["B"] * dis0_cod$B +
efectos_reales_23_vec["C"] * dis0_cod$C +
efectos_reales_23_vec["AB"] * dis0_cod$A * dis0_cod$B +
efectos_reales_23_vec["AC"] * dis0_cod$A * dis0_cod$C +
efectos_reales_23_vec["BC"] * dis0_cod$B * dis0_cod$C +
efectos_reales_23_vec["ABC"] * dis0_cod$A * dis0_cod$B * dis0_cod$C +
error
## Asignación de la respuesta al dis0_cod
## (Diseño factorial completo con réplicas y con 5 puntos centrales)
dis0_cod <- cbind(dis0_cod, yrta)
## Diseño fraccional con réplicas y puntos centrales
dis1_cod <- dis0_cod[c(2,3,5,8,10,11,13,16,17:21),]
## Diseño fraccional con réplicas y sin puntos centrales
dis2_cod <- dis0_cod[c(2,3,5,8,10,11,13,16),]
## Diseño fraccional sin réplica y con cinco puntos centrales
rep1 <- c(9:16,1,4,6,7) #filas que deben ser eliminadas en "sel"
rep2 <- c(1:8,9,12,14,15) #filas que deben ser eliminadas en "sel"
seleccion_name <- sample(c("rep1","rep2"),1) #repetición a eliminar
sel <- get(seleccion_name)
dis3_cod <- dis0_cod[-sel,]
## Diseño fraccional sin réplica y con dos puntos centrales
rep1 <- c(9:16,1,4,6,7) #filas que deben ser eliminadas en "sel"
rep2 <- c(1:8,9,12,14,15) #filas que deben ser eliminadas en "sel"
seleccion_name <- sample(c("rep1","rep2"),1) #repetición a eliminar
sel <- get(seleccion_name)
rep_pc <- sample(c(17:21),3) #punto central a eliminar
dis4_cod <- dis0_cod[-c(sel,rep_pc),]
# Modelo
## Evaluación del diseño con réplicas y con cinco puntos centrales
mod1 <- lm(yrta~A+B+C, data=dis1_cod)
## Evaluación del diseño con réplicas y sin puntos centrales
mod2 <- lm(yrta~A+B+C, data=dis2_cod)
## Evaluación del diseño sin réplicas y con cinco puntos centrales
mod3 <- lm(yrta~A+B+C, data=dis3_cod)
## Evaluación del diseño sin réplicas y con dos puntos centrales
mod4 <- lm(yrta~A+B+C, data=dis4_cod)
## Almacenar los ECM
ECM_tab_23$dis1[i] <- (deviance(mod1) / df.residual(mod1))
ECM_tab_23$dis2[i] <- (deviance(mod2) / df.residual(mod2))
ECM_tab_23$dis3[i] <- (deviance(mod3) / df.residual(mod3))
ECM_tab_23$dis4[i] <- (deviance(mod4) / df.residual(mod4))
# Lista de modelos y de tablas
mods <- list(mod1, mod2, mod3, mod4)
tabs <- list(TAB1_23, TAB2_23, TAB3_23, TAB4_23)
for (j in seq_along(mods)) {
tabs[[j]][i, ] <- sapply(1:3, function(k) {
as.integer(anova(mods[[j]])[k, 5] <= 0.05) })}
## Evaluación del diseño con réplicas y con cinco puntos centrales
TAB1_23 <- tabs[[1]]
## Evaluación del diseño con réplicas y sin puntos centrales
TAB2_23 <- tabs[[2]]
## Evaluación del diseño sin réplicas y con cinco puntos centrales
TAB3_23 <- tabs[[3]]
## Evaluación del diseño sin réplicas y con dos puntos centrales
TAB4_23 <- tabs[[4]]} 4.1.2 Potencia y error tipo I
Código
## Evaluación del diseño con réplicas y con cinco puntos centrales
comp1_23 <- round(colSums(efectos_reales_23[,1:3] == TAB1_23) / n * 100, 2)
pot1_23 <- round(colSums((efectos_reales_23[,1:3] == 1) & (TAB1_23 == 1)) /
colSums(efectos_reales_23[,1:3] == 1) * 100, 2)
err1_23 <- round(colSums((efectos_reales_23[,1:3] == 0) & (TAB1_23 == 1)) /
colSums(efectos_reales_23[,1:3] == 0) * 100, 2)
## Evaluación del diseño con réplicas y sin puntos centrales
comp2_23 <- round(colSums(efectos_reales_23[,1:3] == TAB2_23) / n * 100, 2)
pot2_23 <- round(colSums((efectos_reales_23[,1:3] == 1) & (TAB2_23 == 1)) /
colSums(efectos_reales_23[,1:3] == 1) * 100, 2)
err2_23 <- round(colSums((efectos_reales_23[,1:3] == 0) & (TAB2_23 == 1)) /
colSums(efectos_reales_23[,1:3] == 0) * 100, 2)
## Evaluación del diseño sin réplicas y con cinco puntos centrales
comp3_23 <- round(colSums(efectos_reales_23[,1:3] == TAB3_23) / n * 100, 2)
pot3_23 <- round(colSums((efectos_reales_23[,1:3] == 1) & (TAB3_23 == 1)) /
colSums(efectos_reales_23[,1:3] == 1) * 100, 2)
err3_23 <- round(colSums((efectos_reales_23[,1:3] == 0) & (TAB3_23 == 1)) /
colSums(efectos_reales_23[,1:3] == 0) * 100, 2)
## Evaluación del diseño sin réplicas y con dos puntos centrales
comp4_23 <- round(colSums(efectos_reales_23[,1:3] == TAB4_23) / n * 100, 2)
pot4_23 <- round(colSums((efectos_reales_23[,1:3] == 1) & (TAB4_23 == 1)) /
colSums(efectos_reales_23[,1:3] == 1) * 100, 2)
err4_23 <- round(colSums((efectos_reales_23[,1:3] == 0) & (TAB4_23 == 1)) /
colSums(efectos_reales_23[,1:3] == 0) * 100, 2)4.2 Diseño \(2^{4-1}\)
4.2.1 Significancia de los efectos
Código
# Matrices
TAB1_24 <- matrix(0, nrow = n, ncol = 4)
TAB2_24 <- matrix(0, nrow = n, ncol = 4)
TAB3_24 <- matrix(0, nrow = n, ncol = 4)
TAB4_24 <- matrix(0, nrow = n, ncol = 4)
efectos_reales_24 <- matrix(0, nrow = n, ncol = 15)
names_24 <- c("Ap", "Bp", "Cp", "Dp") # Los que serán estimados por el modelo
colnames(TAB1_24) <- names_24
colnames(TAB2_24) <- names_24
colnames(TAB3_24) <- names_24
colnames(TAB4_24) <- names_24
colnames(efectos_reales_24) <- c("A", "B", "C", "D", #Todos
"AB", "AC", "BC", "AD", "BD", "CD",
"ABC", "ABD", "ACD", "BCD",
"ABCD")
# Tabla para completar los errores
ECM_tab_24 <- data.frame(
dis1 = numeric(n),
dis2 = numeric(n),
dis3 = numeric(n),
dis4 = numeric(n))
# Bucle principal
set.seed(3)
for(i in seq_len(n)){
## Efectos a evaluar en el este diseño
todos_efectos_24 <- c("A", "B", "C", "D",
"AB", "AC", "BC", "AD", "BD", "CD",
"ABC", "ABD", "ACD", "BCD",
"ABCD")
# Selección de un número aleatorio entre 1 y 15
# dado que existen 15 efectos posibles a evaluar en este diseño.
n_ef <- sample(1:length(todos_efectos_24), 1)
# Elección aleatoria _n_ef_ efectos que serán los efectos significativos
efectos_activos <- sample(todos_efectos_24, n_ef)
efectos_reales_24_vec <- setNames(rep(0, length(todos_efectos_24)), todos_efectos_24)
# Asignación de un valor aleatorio a los efectos elegidos como significativos:
efectos_reales_24_vec[efectos_activos] <- runif(n_ef, -7, 7)
# Guardar efectos reales como 1/0
efectos_reales_24[i, ] <- c(
as.integer(efectos_reales_24_vec["A"] != 0),
as.integer(efectos_reales_24_vec["B"] != 0),
as.integer(efectos_reales_24_vec["C"] != 0),
as.integer(efectos_reales_24_vec["D"] != 0),
as.integer(efectos_reales_24_vec["AB"] != 0),
as.integer(efectos_reales_24_vec["AC"] != 0),
as.integer(efectos_reales_24_vec["BC"] != 0),
as.integer(efectos_reales_24_vec["AD"] != 0),
as.integer(efectos_reales_24_vec["BD"] != 0),
as.integer(efectos_reales_24_vec["CD"] != 0),
as.integer(efectos_reales_24_vec["ABC"] != 0),
as.integer(efectos_reales_24_vec["ABD"] != 0),
as.integer(efectos_reales_24_vec["ACD"] != 0),
as.integer(efectos_reales_24_vec["BCD"] != 0),
as.integer(efectos_reales_24_vec["ABCD"] != 0))
# Generar respuesta
## Diseño factorial completo con réplicas y con 5 puntos centrales
dis0_cod <- data.frame(
A = c(rep(c(-1,1), 8), rep(c(-1,1), 8), rep(0,5)),
B = c(rep(c(rep(-1,2), rep(1,2)), 4),
rep(c(rep(-1,2), rep(1,2)), 4),
rep(0,5)),
C = c(rep(c(rep(-1,4), rep(1,4)), 2),
rep(c(rep(-1,4), rep(1,4)), 2),
rep(0,5)),
D = c(rep(-1,8), rep(1,8),
rep(-1,8), rep(1,8),
rep(0,5)))
n_corridas <- nrow(dis0_cod)
## Generación del error
error <- rnorm(n_corridas, 0, 2) #error ~ (mu = 0, sigma^2 = 2^2)
## Generación de la respuesta
yrta <- 50 +
efectos_reales_24_vec["A"] * dis0_cod$A +
efectos_reales_24_vec["B"] * dis0_cod$B +
efectos_reales_24_vec["C"] * dis0_cod$C +
efectos_reales_24_vec["D"] * dis0_cod$D +
efectos_reales_24_vec["AB"] * dis0_cod$A * dis0_cod$B +
efectos_reales_24_vec["AC"] * dis0_cod$A * dis0_cod$C +
efectos_reales_24_vec["BC"] * dis0_cod$B * dis0_cod$C +
efectos_reales_24_vec["AD"] * dis0_cod$A * dis0_cod$D +
efectos_reales_24_vec["BD"] * dis0_cod$B * dis0_cod$D +
efectos_reales_24_vec["CD"] * dis0_cod$C * dis0_cod$D +
efectos_reales_24_vec["ABC"] * dis0_cod$A * dis0_cod$B * dis0_cod$C +
efectos_reales_24_vec["ABD"] * dis0_cod$A * dis0_cod$B * dis0_cod$D +
efectos_reales_24_vec["ACD"] * dis0_cod$A * dis0_cod$C * dis0_cod$D +
efectos_reales_24_vec["BCD"] * dis0_cod$B * dis0_cod$C * dis0_cod$D +
efectos_reales_24_vec["ABCD"] * dis0_cod$A * dis0_cod$B * dis0_cod$C * dis0_cod$D +
error
## Asignación de la respuesta al dis0_cod
## (Diseño factorial completo con réplicas y con 5 puntos centrales)
dis0_cod <- cbind(dis0_cod, yrta)
## Diseño fraccional con réplicas y puntos centrales
dis1_cod <- dis0_cod[c(1,4,6,7,10,11,13,16,
17,20,22,23,26,27,29,32,
33:37),]
## Diseño con réplicas y sin puntos centrales
dis2_cod <- dis0_cod[c(1,4,6,7,10,11,13,16,
17,20,22,23,26,27,29,32),]
## Diseño sin réplica y con cinco puntos centrales
rep1 <- c(17:32,2,3,5,8,9,12,14,15) #filas para eliminar en "sel"
rep2 <- c(1:16,18,19,21,24,25,28,30,31) #filas para eliminar en "sel"
seleccion_name <- sample(c("rep1","rep2"),1) #repetición a eliminar
sel <- get(seleccion_name)
dis3_cod <- dis0_cod[-sel,]
## Diseño sin réplica y con dos puntos centrales
rep1 <- c(17:32,2,3,5,8,9,12,14,15) #filas para eliminar en "sel"
rep2 <- c(1:16,18,19,21,24,25,28,30,31) #filas para eliminar en "sel"
seleccion_name <- sample(c("rep1","rep2"),1) #repetición a eliminar
sel <- get(seleccion_name)
rep_pc <- sample(c(33:37),3) #punto central a eliminar
dis4_cod <- dis0_cod[-c(sel,rep_pc),]
# Modelo
## Evaluación del diseño con réplicas y con cinco puntos centrales
mod1 <- lm(yrta ~ A + B + C + D, data=dis1_cod)
## Evaluación del diseño con réplicas y sin puntos centrales
mod2 <- lm(yrta ~ A + B + C + D, data=dis2_cod)
## Evaluación del diseño sin réplicas y con cinco puntos centrales
mod3 <- lm(yrta ~ A + B + C + D, data=dis3_cod)
## Evaluación del diseño sin réplicas y con dos puntos centrales
mod4 <- lm(yrta ~ A + B + C + D, data=dis4_cod)
## Almacenar los ECM
ECM_tab_24$dis1[i] <- (deviance(mod1) / df.residual(mod1))
ECM_tab_24$dis2[i] <- (deviance(mod2) / df.residual(mod2))
ECM_tab_24$dis3[i] <- (deviance(mod3) / df.residual(mod3))
ECM_tab_24$dis4[i] <- (deviance(mod4) / df.residual(mod4))
# Lista de modelos y de tablas
mods <- list(mod1, mod2, mod3, mod4)
tabs <- list(TAB1_24, TAB2_24, TAB3_24, TAB4_24)
for (j in seq_along(mods)) {
tabs[[j]][i, ] <- sapply(1:4, function(k) {
as.integer(anova(mods[[j]])[k, 5] <= 0.05) })}
## Diseño con réplicas y con cinco puntos centrales
TAB1_24 <- tabs[[1]]
## Evaluación del diseño con réplicas y sin puntos centrales
TAB2_24 <- tabs[[2]]
## Evaluación del diseño sin réplicas y con cinco puntos centrales
TAB3_24 <- tabs[[3]]
## Evaluación del diseño sin réplicas y con dos puntos centrales
TAB4_24 <- tabs[[4]]} 4.2.2 Potencia y error tipo I
Código
## Evaluación del diseño con réplicas y con cinco puntos centrales
comp1_24 <- round(colSums(efectos_reales_24[,1:4] == TAB1_24) / n * 100, 2)
pot1_24 <- round(colSums((efectos_reales_24[,1:4] == 1) & (TAB1_24 == 1)) /
colSums(efectos_reales_24[,1:4] == 1) * 100, 2)
err1_24 <- round(colSums((efectos_reales_24[,1:4] == 0) & (TAB1_24 == 1)) /
colSums(efectos_reales_24[,1:4] == 0) * 100, 2)
## Evaluación del diseño con réplicas y sin puntos centrales
comp2_24 <- round(colSums(efectos_reales_24[,1:4] == TAB2_24) / n * 100, 2)
pot2_24 <- round(colSums((efectos_reales_24[,1:4] == 1) & (TAB2_24 == 1)) /
colSums(efectos_reales_24[,1:4] == 1) * 100, 2)
err2_24 <- round(colSums((efectos_reales_24[,1:4] == 0) & (TAB2_24 == 1)) /
colSums(efectos_reales_24[,1:4] == 0) * 100, 2)
## Evaluación del diseño sin réplicas y con cinco puntos centrales
comp3_24 <- round(colSums(efectos_reales_24[,1:4] == TAB3_24) / n * 100, 2)
pot3_24 <- round(colSums((efectos_reales_24[,1:4] == 1) & (TAB3_24 == 1)) /
colSums(efectos_reales_24[,1:4] == 1) * 100, 2)
err3_24 <- round(colSums((efectos_reales_24[,1:4] == 0) & (TAB3_24 == 1)) /
colSums(efectos_reales_24[,1:4] == 0) * 100, 2)
## Evaluación del diseño sin réplicas y con dos puntos centrales
comp4_24 <- round(colSums(efectos_reales_24[,1:4] == TAB4_24) / n * 100, 2)
pot4_24 <- round(colSums((efectos_reales_24[,1:4] == 1) & (TAB4_24 == 1)) /
colSums(efectos_reales_24[,1:4] == 1) * 100, 2)
err4_24 <- round(colSums((efectos_reales_24[,1:4] == 0) & (TAB4_24 == 1)) /
colSums(efectos_reales_24[,1:4] == 0) * 100, 2)4.3 Diseño \(2^{5-2}\)
4.3.1 Significancia de los efectos
Código
# Matrices
TAB1_25 <- matrix(0, nrow = n, ncol = 5)
TAB2_25 <- matrix(0, nrow = n, ncol = 5)
TAB3_25 <- matrix(0, nrow = n, ncol = 5)
TAB4_25 <- matrix(0, nrow = n, ncol = 5)
efectos_reales_25 <- matrix(0, nrow = n, ncol = 31)
names_25 <- c("Ap","Bp","Cp","Dp","Ep")
colnames(TAB1_25) <- names_25
colnames(TAB2_25) <- names_25
colnames(TAB3_25) <- names_25
colnames(TAB4_25) <- names_25
colnames(efectos_reales_25) <- c("A","B","C","D","E",
"AB","AC","BC","AD","BD","CD","AE","BE","CE","DE",
"ABC","ABD","ACD","BCD","ABE","ACE","BCE","ADE","BDE","CDE",
"ABCD","ABCE","ABDE","ACDE","BCDE",
"ABCDE")
# Tabla para completar los errores
ECM_tab_25 <- data.frame(
dis1 = numeric(n),
dis2 = numeric(n),
dis3 = numeric(n),
dis4 = numeric(n))
# Bucle principal
set.seed(3)
for(i in seq_len(n)){
## Efectos a evaluar en el este diseño
todos_efectos_25 <- c("A","B","C","D","E",
"AB","AC","BC","AD","BD","CD","AE","BE","CE","DE",
"ABC","ABD","ACD","BCD","ABE","ACE","BCE","ADE","BDE","CDE",
"ABCD","ABCE","ABDE","ACDE","BCDE",
"ABCDE")
# Selección de un número aleatorio entre 1 y 31
# dado que existen 31 efectos posibles a evaluar en este diseño.
n_ef <- sample(1:length(todos_efectos_25), 1)
# Elección aleatoria _n_ef_ efectos que serán los efectos significativos
efectos_activos <- sample(todos_efectos_25, n_ef)
efectos_reales_25_vec <- setNames(rep(0, length(todos_efectos_25)), todos_efectos_25)
# Asignación de un valor aleatorio a los efectos elegidos como significativos:
efectos_reales_25_vec[efectos_activos] <- sample(c(runif(n_ef, -7, -0.5),
runif(n_ef, 0.5, 7)), n_ef)
# Guardar efectos reales como 1/0
efectos_reales_25[i, ] <- c(
as.integer(efectos_reales_25_vec["A"] != 0),
as.integer(efectos_reales_25_vec["B"] != 0),
as.integer(efectos_reales_25_vec["C"] != 0),
as.integer(efectos_reales_25_vec["D"] != 0),
as.integer(efectos_reales_25_vec["E"] != 0),
as.integer(efectos_reales_25_vec["AB"] != 0),
as.integer(efectos_reales_25_vec["AC"] != 0),
as.integer(efectos_reales_25_vec["BC"] != 0),
as.integer(efectos_reales_25_vec["AD"] != 0),
as.integer(efectos_reales_25_vec["BD"] != 0),
as.integer(efectos_reales_25_vec["CD"] != 0),
as.integer(efectos_reales_25_vec["AE"] != 0),
as.integer(efectos_reales_25_vec["BE"] != 0),
as.integer(efectos_reales_25_vec["CE"] != 0),
as.integer(efectos_reales_25_vec["DE"] != 0),
as.integer(efectos_reales_25_vec["ABC"] != 0),
as.integer(efectos_reales_25_vec["ABD"] != 0),
as.integer(efectos_reales_25_vec["ACD"] != 0),
as.integer(efectos_reales_25_vec["BCD"] != 0),
as.integer(efectos_reales_25_vec["ABE"] != 0),
as.integer(efectos_reales_25_vec["ACE"] != 0),
as.integer(efectos_reales_25_vec["BCE"] != 0),
as.integer(efectos_reales_25_vec["ADE"] != 0),
as.integer(efectos_reales_25_vec["BDE"] != 0),
as.integer(efectos_reales_25_vec["CDE"] != 0),
as.integer(efectos_reales_25_vec["ABCD"] != 0),
as.integer(efectos_reales_25_vec["ABCE"] != 0),
as.integer(efectos_reales_25_vec["ABDE"] != 0),
as.integer(efectos_reales_25_vec["ACDE"] != 0),
as.integer(efectos_reales_25_vec["BCDE"] != 0),
as.integer(efectos_reales_25_vec["ABCDE"]!= 0))
# Generar respuesta
## Diseño de variables codificadas con réplicas y con 5 puntos centrales
dis0_cod <- data.frame(
A = c(rep(c(-1,1),16), rep(c(-1,1),16), rep(0,5)),
B = c(rep(c(rep(-1,2), rep(1,2)),8), rep(c(rep(-1,2), rep(1,2)),8),rep(0,5)),
C = c(rep(c(rep(-1,4), rep(1,4)),4), rep(c(rep(-1,4), rep(1,4)),4),rep(0,5)),
D = c(rep(c(rep(-1,8), rep(1,8)),2), rep(c(rep(-1,8), rep(1,8)),2),rep(0,5)),
E = c(rep(-1,16), rep(1,16), rep(-1,16), rep(1,16), rep(0,5)))
n_corridas <- nrow(dis0_cod)
## Generación del error
error <- rnorm(n_corridas, 0, 2) #error ~ (mu = 0, sigma^2 = 2^2)
## Generación de la respuesta
yrta <- 50 +
efectos_reales_25_vec["A"] * dis0_cod$A +
efectos_reales_25_vec["B"] * dis0_cod$B +
efectos_reales_25_vec["C"] * dis0_cod$C +
efectos_reales_25_vec["D"] * dis0_cod$D +
efectos_reales_25_vec["E"] * dis0_cod$E +
efectos_reales_25_vec["AB"] * dis0_cod$A * dis0_cod$B +
efectos_reales_25_vec["AC"] * dis0_cod$A * dis0_cod$C +
efectos_reales_25_vec["BC"] * dis0_cod$B * dis0_cod$C +
efectos_reales_25_vec["AD"] * dis0_cod$A * dis0_cod$D +
efectos_reales_25_vec["BD"] * dis0_cod$B * dis0_cod$D +
efectos_reales_25_vec["CD"] * dis0_cod$C * dis0_cod$D +
efectos_reales_25_vec["AE"] * dis0_cod$A * dis0_cod$E +
efectos_reales_25_vec["BE"] * dis0_cod$B * dis0_cod$E +
efectos_reales_25_vec["CE"] * dis0_cod$C * dis0_cod$E +
efectos_reales_25_vec["DE"] * dis0_cod$D * dis0_cod$E +
efectos_reales_25_vec["ABC"] * dis0_cod$A * dis0_cod$B * dis0_cod$C +
efectos_reales_25_vec["ABD"] * dis0_cod$A * dis0_cod$B * dis0_cod$D +
efectos_reales_25_vec["ACD"] * dis0_cod$A * dis0_cod$C * dis0_cod$D +
efectos_reales_25_vec["BCD"] * dis0_cod$B * dis0_cod$C * dis0_cod$D +
efectos_reales_25_vec["ABE"] * dis0_cod$A * dis0_cod$B * dis0_cod$E +
efectos_reales_25_vec["ACE"] * dis0_cod$A * dis0_cod$C * dis0_cod$E +
efectos_reales_25_vec["BCE"] * dis0_cod$B * dis0_cod$C * dis0_cod$E +
efectos_reales_25_vec["ADE"] * dis0_cod$A * dis0_cod$D * dis0_cod$E +
efectos_reales_25_vec["BDE"] * dis0_cod$B * dis0_cod$D * dis0_cod$E +
efectos_reales_25_vec["CDE"] * dis0_cod$C * dis0_cod$D * dis0_cod$E +
efectos_reales_25_vec["ABCD"] * dis0_cod$A * dis0_cod$B * dis0_cod$C * dis0_cod$D +
efectos_reales_25_vec["ABCE"] * dis0_cod$A * dis0_cod$B * dis0_cod$C * dis0_cod$E +
efectos_reales_25_vec["ABDE"] * dis0_cod$A * dis0_cod$B * dis0_cod$D * dis0_cod$E +
efectos_reales_25_vec["ACDE"] * dis0_cod$A * dis0_cod$C * dis0_cod$D * dis0_cod$E +
efectos_reales_25_vec["BCDE"] * dis0_cod$B * dis0_cod$C * dis0_cod$D * dis0_cod$E +
efectos_reales_25_vec["ABCDE"] * dis0_cod$A * dis0_cod$B * dis0_cod$C * dis0_cod$D * dis0_cod$E +
error
## Asignación de la respuesta al dis0_cod
dis0_cod <- cbind(dis0_cod, yrta)
## Diseño fraccional con réplicas y puntos centrales
dis1_cod <- dis0_cod[c(2, 3, 5, 8, 9, 12, 14, 15, 17, 20, 22, 23, 26, 27, 29, 32,
34, 35, 37, 40, 41, 44, 46, 47, 49, 52, 54, 55, 58, 59, 61, 64,
65:69),]
## Diseño con réplicas y sin puntos centrales
dis2_cod <- dis0_cod[c(2, 3, 5, 8, 9, 12, 14, 15, 17, 20, 22, 23, 26, 27, 29, 32,
34, 35, 37, 40, 41, 44, 46, 47, 49, 52, 54, 55, 58, 59, 61, 64),]
## Diseño sin réplica y con cinco puntos centrales
rep1 <- c(33:64, 1, 4, 6, 7, 10, 11, 13, 16, 18, 19, 21, 24, 25, 28, 30, 31)
rep2 <- c(1:32, 33, 36, 38, 39, 42, 43, 45, 48, 50, 51, 53, 56, 57, 60, 62, 63)
seleccion_name <- sample(c("rep1","rep2"),1) #repetición a eliminar
sel <- get(seleccion_name)
dis3_cod <- dis0_cod[-sel,]
## Diseño sin réplica y con dos puntos centrales
rep1 <- c(33:64, 1, 4, 6, 7, 10, 11, 13, 16, 18, 19, 21, 24, 25, 28, 30, 31)
rep2 <- c(1:32, 33, 36, 38, 39, 42, 43, 45, 48, 50, 51, 53, 56, 57, 60, 62, 63)
seleccion_name <- sample(c("rep1","rep2"),1) #repetición a eliminar
sel <- get(seleccion_name)
rep_pc <- sample(c(65:69),3) #punto central a eliminar
dis4_cod <- dis0_cod[-c(sel,rep_pc),]
# Modelo
## Evaluación del diseño con réplicas y con cinco puntos centrales
mod1 <- lm(yrta ~ A + B + C + D + E, data=dis1_cod)
## Evaluación del diseño con réplicas y sin puntos centrales
mod2 <- lm(yrta ~ A + B + C + D + E, data=dis2_cod)
## Evaluación del diseño sin réplicas y con cinco puntos centrales
mod3 <- lm(yrta ~ A + B + C + D + E, data=dis3_cod)
## Evaluación del diseño sin réplicas y con dos puntos centrales
mod4 <- lm(yrta ~ A + B + C + D + E, data=dis4_cod)
## Almacenar los ECM
ECM_tab_25$dis1[i] <- (deviance(mod1) / df.residual(mod1))
ECM_tab_25$dis2[i] <- (deviance(mod2) / df.residual(mod2))
ECM_tab_25$dis3[i] <- (deviance(mod3) / df.residual(mod3))
ECM_tab_25$dis4[i] <- (deviance(mod4) / df.residual(mod4))
# Lista de modelos y de tablas
mods <- list(mod1, mod2, mod3, mod4)
tabs <- list(TAB1_25, TAB2_25, TAB3_25, TAB4_25)
for (j in seq_along(mods)) {
tabs[[j]][i, ] <- sapply(1:5, function(k) {
as.integer(anova(mods[[j]])[k, 5] <= 0.05) })}
## Diseño con réplicas y con cinco puntos centrales
TAB1_25 <- tabs[[1]]
## Evaluación del diseño con réplicas y sin puntos centrales
TAB2_25 <- tabs[[2]]
## Evaluación del diseño sin réplicas y con cinco puntos centrales
TAB3_25 <- tabs[[3]]
## Evaluación del diseño sin réplicas y con dos puntos centrales
TAB4_25 <- tabs[[4]]} 4.3.2 Potencia y error tipo I
Código
## Evaluación del diseño con réplicas y con cinco puntos centrales
comp1_25 <- round(colSums(efectos_reales_25[,c(1:5)] == TAB1_25) / n * 100, 2)
pot1_25 <- round(colSums((efectos_reales_25[,c(1:5)] == 1) & (TAB1_25 == 1)) /
colSums(efectos_reales_25[,c(1:5)] == 1) * 100, 2)
err1_25 <- round(colSums((efectos_reales_25[,c(1:5)] == 0) & (TAB1_25 == 1)) /
colSums(efectos_reales_25[,c(1:5)] == 0) * 100, 2)
## Evaluación del diseño con réplicas y sin puntos centrales
comp2_25 <- round(colSums(efectos_reales_25[,c(1:5)] == TAB2_25) / n * 100, 2)
pot2_25 <- round(colSums((efectos_reales_25[,c(1:5)] == 1) & (TAB2_25 == 1)) /
colSums(efectos_reales_25[,c(1:5)] == 1) * 100, 2)
err2_25 <- round(colSums((efectos_reales_25[,c(1:5)] == 0) & (TAB2_25 == 1)) /
colSums(efectos_reales_25[,c(1:5)] == 0) * 100, 2)
## Evaluación del diseño sin réplicas y con cinco puntos centrales
comp3_25 <- round(colSums(efectos_reales_25[,c(1:5)] == TAB3_25) / n * 100, 2)
pot3_25 <- round(colSums((efectos_reales_25[,c(1:5)] == 1) & (TAB3_25 == 1)) /
colSums(efectos_reales_25[,c(1:5)] == 1) * 100, 2)
err3_25 <- round(colSums((efectos_reales_25[,c(1:5)] == 0) & (TAB3_25 == 1)) /
colSums(efectos_reales_25[,c(1:5)] == 0) * 100, 2)
## Evaluación del diseño sin réplicas y con dos puntos centrales
comp4_25 <- round(colSums(efectos_reales_25[,c(1:5)] == TAB4_25) / n * 100, 2)
pot4_25 <- round(colSums((efectos_reales_25[,c(1:5)] == 1) & (TAB4_25 == 1)) /
colSums(efectos_reales_25[,c(1:5)] == 1) * 100, 2)
err4_25 <- round(colSums((efectos_reales_25[,c(1:5)] == 0) & (TAB4_25 == 1)) /
colSums(efectos_reales_25[,c(1:5)] == 0) * 100, 2)4.4 Evaluación del Error Cuadrático Medio (ECM) por medidas repetidas
Para la evaluación del Error Cuadrático Medio (ECM) se diseñó la siguiente función:
Código
# X: Matriz que incluye los datos para comparar que debe estar organizada así
# C: Matriz de contrastes
# alpha: Nivel de significancia
Hotelling <- function(X, C, alpha) {
## Cálculo de la $T^2$ Hotelling
mu <- colMeans(X)
n <- nrow(X) # Número de individuos (filas de X)
p <- ncol(X) # Número de tratamientos (columnas de X)
# Matriz identidad e I - (1/n)11'
In <- diag(n)
uno <- matrix(1, n, n) # 1 x t(1)
H <- In - (1/n) * uno
# Calculo de S
S <- (1/(n-1)) * t(X) %*% H %*% X
# Cálculo de T2
T2 <- n*(t(C%*%mu))%*%(solve(C%*%S%*%t(C)))%*%(C%*%mu)
# F calculado
Fcal <- (n - p + 1) / ((n - 1)*(p - 1)) * T2
# F tabulado (crítico)
gl1 <- p-1
gl2 <- n-p+1
# Valor crítico F(0.05; p-1, n-p+1)
F_critico <- qf(1-alpha, df1 = gl1, df2 = gl2)
return(list(
T2 = as.numeric(T2),
Fcal = as.numeric(Fcal),
F_critico = F_critico))
}5 Resultados
5.1 Diseño \(2^{3-1}\)
Factores | Con réplicas con 5 puntos centrales | Con réplicas con 5 puntos centrales | Con réplicas con 5 puntos centrales | Con réplicas sin puntos centrales | Con réplicas sin puntos centrales | Con réplicas sin puntos centrales | Sin réplicas con 5 puntos centrales | Sin réplicas con 5 puntos centrales | Sin réplicas con 5 puntos centrales | Sin réplicas con 2 puntos centrales | Sin réplicas con 2 puntos centrales | Sin réplicas con 2 puntos centrales |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
Concordancia | Potencia | Error Tipo I | Concordancia | Potencia | Error Tipo I | Concordancia | Potencia | Error Tipo I | Concordancia | Potencia | Error Tipo I | |
Ap | 50 | 40.00 | 40.00 | 50 | 40.00 | 40.00 | 50 | 40.00 | 40.00 | 60 | 40.00 | 20 |
Bp | 80 | 80.00 | 20.00 | 80 | 80.00 | 20.00 | 70 | 60.00 | 20.00 | 70 | 40.00 | 0 |
Cp | 70 | 66.67 | 28.57 | 60 | 33.33 | 28.57 | 60 | 33.33 | 28.57 | 80 | 33.33 | 0 |
5.2 Diseño \(2^{4-1}\)
Factores | Con réplicas con 5 puntos centrales | Con réplicas con 5 puntos centrales | Con réplicas con 5 puntos centrales | Con réplicas sin puntos centrales | Con réplicas sin puntos centrales | Con réplicas sin puntos centrales | Sin réplicas con 5 puntos centrales | Sin réplicas con 5 puntos centrales | Sin réplicas con 5 puntos centrales | Sin réplicas con 2 puntos centrales | Sin réplicas con 2 puntos centrales | Sin réplicas con 2 puntos centrales |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
Concordancia | Potencia | Error Tipo I | Concordancia | Potencia | Error Tipo I | Concordancia | Potencia | Error Tipo I | Concordancia | Potencia | Error Tipo I | |
Ap | 40 | 0.00 | 20.0 | 40 | 0.00 | 20.0 | 40 | 0 | 20.0 | 50 | 0 | 0.0 |
Bp | 50 | 33.33 | 25.0 | 40 | 16.67 | 25.0 | 30 | 0 | 25.0 | 40 | 0 | 0.0 |
Cp | 80 | 50.00 | 12.5 | 80 | 50.00 | 12.5 | 70 | 0 | 12.5 | 70 | 0 | 12.5 |
Dp | 40 | 20.00 | 40.0 | 40 | 20.00 | 40.0 | 50 | 20 | 20.0 | 50 | 20 | 20.0 |
5.3 Diseño \(2^5\)
Factores | Con réplicas con 5 puntos centrales | Con réplicas con 5 puntos centrales | Con réplicas con 5 puntos centrales | Con réplicas sin puntos centrales | Con réplicas sin puntos centrales | Con réplicas sin puntos centrales | Sin réplicas con 5 puntos centrales | Sin réplicas con 5 puntos centrales | Sin réplicas con 5 puntos centrales | Sin réplicas con 2 puntos centrales | Sin réplicas con 2 puntos centrales | Sin réplicas con 2 puntos centrales |
|---|---|---|---|---|---|---|---|---|---|---|---|---|
Concordancia | Potencia | Error Tipo I | Concordancia | Potencia | Error Tipo I | Concordancia | Potencia | Error Tipo I | Concordancia | Potencia | Error Tipo I | |
Ap | 70 | 50.00 | 16.67 | 80 | 50.00 | 0.00 | 60 | 0.00 | 0.00 | 60 | 0.00 | 0.00 |
Bp | 40 | 0.00 | 33.33 | 40 | 0.00 | 33.33 | 40 | 0.00 | 33.33 | 50 | 0.00 | 16.67 |
Cp | 60 | 25.00 | 16.67 | 60 | 25.00 | 16.67 | 60 | 0.00 | 0.00 | 60 | 0.00 | 0.00 |
Dp | 80 | 66.67 | 14.29 | 80 | 66.67 | 14.29 | 90 | 66.67 | 0.00 | 90 | 66.67 | 0.00 |
Ep | 70 | 60.00 | 20.00 | 70 | 60.00 | 20.00 | 70 | 60.00 | 20.00 | 50 | 20.00 | 20.00 |
5.4 Evaluación del ECM por medidas repetidas para todos los tratamientos
5.4.1 Análisis descriptivo de los ECM por diseño
5.4.2 Matriz de contrastes generales \(C_1\)
Diseño | T² | F_calculado | F_crítico | Decisión |
|---|---|---|---|---|
Diseño factorial 2³⁻¹ | 1.34 | 0.35 | 4.35 | No se rechaza H₀ |
Diseño factorial 2⁴⁻¹ | 20.53 | 5.32 | 4.35 | Se rechaza H₀ |
Diseño factorial 2⁵⁻² | 18.62 | 4.83 | 4.35 | Se rechaza H₀ |
5.4.3 Matriz de contrastes generales de interés \(C_2\)
Diseño | T² | F_calculado | F_crítico | Decisión |
|---|---|---|---|---|
Diseño factorial 2³⁻¹ | 32.93 | 8.54 | 4.35 | Se rechaza H₀ |
Diseño factorial 2⁴⁻¹ | 21.59 | 5.60 | 4.35 | Se rechaza H₀ |
Diseño factorial 2⁵⁻² | 20.63 | 5.35 | 4.35 | Se rechaza H₀ |
Dado que, al evaluar la matriz de contrastes, se rechazó la hipótesis nula, es decir, al menos hay uno de los contrastes que es diferente de 0, resulta necesario analizar cada contraste de manera individual, obteniéndose los siguientes resultados:
5.4.3.1 Contraste \(C_{21} =\) (-1, -1, 1, 1)
Diseño | T² | F_calculado | F_crítico | Decisión |
|---|---|---|---|---|
Diseño factorial 2³⁻¹ | 0.46 | 0.12 | 4.35 | No se rechaza H₀ |
Diseño factorial 2⁴⁻¹ | 5.87 | 1.52 | 4.35 | No se rechaza H₀ |
Diseño factorial 2⁵⁻² | 7.09 | 1.84 | 4.35 | No se rechaza H₀ |
5.4.3.2 Contraste \(C_{22} =\) (-1, 1, 1, 1)
Diseño | T² | F_calculado | F_crítico | Decisión |
|---|---|---|---|---|
Diseño factorial 2³⁻¹ | 29.53 | 7.66 | 4.35 | Se rechaza H₀ |
Diseño factorial 2⁴⁻¹ | 7.95 | 2.06 | 4.35 | No se rechaza H₀ |
Diseño factorial 2⁵⁻² | 14.02 | 3.63 | 4.35 | No se rechaza H₀ |
5.4.3.3 Contraste \(C_{23} =\) (1, -1, 0, 0)
Diseño | T² | F_calculado | F_crítico | Decisión |
|---|---|---|---|---|
Diseño factorial 2³⁻¹ | 1.05 | 0.27 | 4.35 | No se rechaza H₀ |
Diseño factorial 2⁴⁻¹ | 6.14 | 1.59 | 4.35 | No se rechaza H₀ |
Diseño factorial 2⁵⁻² | 14.28 | 3.70 | 4.35 | No se rechaza H₀ |
Dı́az Monroy, Luis Guillermo, y Rivera Morales Mario Alfonso. 2012. Análisis estadı́stico de datos multivariados. Universidad Nacional de Colombia.